arXiv: 1503.02014vl [cond-mat.soft] 6 Mar 2015 


Ferroelectric domain formation in discotic liquid crystals : Monte Carlo study on the 

influence of boundary conditions 

Tushar Kanti Boscf] and Jayashree SahaQ 

Department of Physics, University of Calcutta, 92 Acharya Prafulla Chandra Road, Kolkata 700009, India 

(Dated: March 9, 2015) 

The realization of a spontaneous macroscopic ferroelectric order in fluids of anisotropic mesogens 
is a topic of both fundamental and technological interest. Recently, we demonstrated that a system 
of dipolar achiral disklike ellipsoids can exhibit long-searched ferroelectric liquid crystalline phases 
of dipolar origin. In the present work, extensive off-lattice Monte Carlo simulations are used to 
investigate the phase behavior of the system under the influences of the electrostatic boundary con¬ 
ditions that restrict any global polarization. We find that the system develops strongly ferroelectric 
slablike domains periodically arranged in an antiferroelectric fashion. Exploring the phase behavior 
at different dipole strengths, we find existence of the ferroelectric nematic and ferroelectric columnar 
order inside the domains. For higher dipole strengths, a biaxial phase is also obtained with a similar 
periodic array of ferroelectric slabs of antiparallel polarizations. We have studied the depolarizing 
effects by using both the Ewald summation and the spherical cut-off techniques. We present and 
compare the results of the two different approaches of considering the depolarizing effects in this 
anisotropic system. It is explicitly shown that the domain size increases with the system size as a 
result of considering longer range of dipolar interactions. The system exhibits pronounced system 
size effects for stronger dipolar interactions. The results provide strong evidence to the novel under¬ 
standing that the dipolar interactions are indeed sufficient to produce long range ferroelectric order 
in anisotropic fluids. 


I. INTRODUCTION 

The understanding and development of various 
possible polar order in liquid crystals is an important 
area of research in soft matter physics and chemistry 
00. The studies of ferroelectric and antiferroelectric 
fluids has attracted much attention in the last 25 years 
[3, |H-@. The systems of chiral rodlike molecules and 
achiral bent core molecules are the conventional systems 
which exhibited ferroelectric and antiferroelectric liquid 
crystal phases. Chiral disk shaped molecules also 
exhibit ferroelectric columnar phases [6]. Usually in the 
conventional ferroelectric smectic and columnar phases, 
the respective ferroelectric orderings are not primarily 
driven by dipole-dipole interactions. However, there is 
no fundamental reason to forbid a proper ferroelectric 
fluid with spontaneous ferroelectric order of dipolar 
origin 0®. From computer simulation studies, it 
has been found that model spherical particles with a 
strong central dipole moment exhibit ferroelectric fluid 
phases in conducting surroundings 0 0 EMI. In 
some studies, the dipolar spheres showed the formation 
of ferroelectric domains in the presence of a strong 
depolarizing field in insulating vacuum surroundings 
OS 17 ]. The spontaneous ferroelectric order in the 
system of dipolar spheres is developed solely due to 
the dipolar interactions and the lack of orientational 
bias of the spherical particles. On the contrary, the 
situation is very different for aspherical particles having 
dipole moments. A spontaneous macroscopic ferro¬ 


electric order of dipolar origin is very rare in fluids of 
anisotropic molecules which are the natural candidates 
to form liquid crystal phases. It was also shown that 
as the length to breadth ratio decreases from unity, the 
tendency to form a ferroelectric nematic phase gradually 
decreases for discotic ellipsoidal particles with strong 
central dipoles 0 - Recently, we have shown that 
the realization of a novel class of proper ferroelectric 
liquid crystal phases is possible in a system of achiral 
dipolar disklike Gay-Berne(GB) ellipsoids in conducting 
surroundings (2(|. The model system of interest consists 
of attractive-repulsive GB oblate ellipsoids embedded 
with two parallel point dipoles positioned symmetrically 
on the equatorial plane of the ellipsoids. The system 
exhibited stable ferroelectric nematic and columnar 
fluids with strong overall polarization. The study 
demonstrated that the dipolar interactions are indeed 
sufficient to produce a class of novel ferroelectric fluids 
of essential interest in systems of disc shaped particles 
poj . A system exhibiting a spontaneous macroscopic 
polarization is usually expected to be sensitive to the 
electrostatic boundary conditions. Thus, it becomes 
interesting enough to investigate the influences of the 
boundary conditions upon the novel ferroelectric liquid 
crystal phases of interest. Here, we study the phase 
behavior of the system under the influence of the bound¬ 
ary condition that restricts any global polarization. We 
study the phase behavior of the system using the simple 
spherical cutoff approach excluding the contribution of 
a polarizing reaction field [2l| and also with the more 
conventional Ewald summation method 0 incorporat¬ 
ing the effect of a strong depolarizing field. Both the 
methods are described in Sec. un 
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Extensive work focused on testing the influence of 
different boundary conditions on the behavior of simple 
models of dipolar liquids, for example, the dipolar hard 
spheres and Stockmayer fluid, was started in the 1980s 
[ 2214271 ] . Even recently computational studies of dipolar 
systems have been performed using different boundary 
conditions mm. However, the present work is one of 
the firsts in testing the effects of boundary conditions 
on an anisotropic system exhibiting spontaneous polar¬ 
ization of dipolar origin. It must be emphasized that 
the attractive-repulsive GB ellipsoids can naturally form 
liquid crystal phases even without dipolar interactions 
whereas the spherically symmetric particles are unable 
to form orientationally ordered phases without dipolar 
interactions. 

Now, it can be easily understood that the problem of 
determining the collective organization of a system of 
dipolar molecules is definitely non-trivial, since it is the 
result of a balance between the tendency of two dipoles to 
arrange anti-parallel and the formation of domains with 
a common dipole orientation that gives fascinating do¬ 
main structures [13, 0, In the present work, we 

systematically investigate the influence of the depolariz¬ 
ing boundary condition upon the existence of different 
ferroelectric phases of interest. The system results in the 
formation of slablike ferroelectric fluid domains periodi¬ 
cally arranged in an antiferroelectric fashion. More im¬ 
portantly, at different state points and dipole strengths, 
we find existence of long range ferroelectric nematic, fer¬ 
roelectric columnar and ferroelectric biaxial order within 
the domains. Systems of different size (N = 1500, 4000 
and 8500) have been investigated via the simple spheri¬ 
cal cutoff approach. It is explicitly shown that the size 
of the ferroelectric domains grow with the system size as 
an effect of considering longer range of dipolar interac¬ 
tions. In Sec. un we detail the molecular model and pair 
interactions. The simulation methods are described in 
Sec. m In Sec. El we describe the simulation results 
and Sec. IVl concludes the paper. 


II. MODEL 

In this paper we present computer simulations of a sys¬ 
tem of uniaxial oblate GB ellipsoids of revolution where 
each ellipsoid is embedded with two axial off-center par¬ 
allel point dipole moments. The dipoles are symmetri¬ 
cally placed on the equatorial plane of the ellipsoid, at 
equal distances from the center of the ellipsoid. In the 
present work, the dipoles are placed on the molecular x 
axis (perpendicular to the symmetry axis) of each GB 
molecule, separated by a distance d* = d/ao = 0.5 along 
the axis. The dipolar ellipsoids are interacting via a pair 
potential which is a sum of a modified form of the GB 
potential 0 and the electrostatic dipolar interactions. 
In the modified form for discotic liquid crystal [32], the 
pair potential between two oblate ellipsoids i and j is 


given by 

uj) = 4e(fy, u i; Uj)(p V 2 - p“ 6 ) 

where pij = (r^ — a(r^-, uq, u j) + a e )/a e . Here unit 
vectors up and u j represent the orientations of the 
symmetry axes of the molecules, r^- = is the 

separation vector of length between the centers of 
mass of the ellipsoids and a e is the minimum separation 
between two ellipsoids in a face-to-face configuration de¬ 
termining the thickness of the ellipsoids. The anisotropic 
contact distance a and the depth of pair interaction 
well e are dependent on four important parameters 
/€, k/, /i, z/, as defined in [32[. Here k = a e /a 0 is the 
aspect ratio of the ellipsoids where ao is the minimum 
separation between two ellipsoids in a side-by-side 
configuration, k' = e s /e e is ratio of interaction well 
depths in side-by-side and face-to-face configuration of 
the disc shaped ellipsoids. The other two parameters 
fi and v control the well depth of the potential. ao 
and eo define the length and energy scales respectively 
where eo is the well depth in the cross configuration. 
The values used here to study the bulk phase behavior 
are k = 0.345, k' = 0.2,/i = 1, z/ = 3. The value of k is 
obtained from the parametrization of the GB potential 
that mimics the interaction between two molecules of 
triphynylene [33] which is known to form the core of 
many discotic mesogens (34[. The other parameters 
were chosen from previous works on discotic liquid 
crystals which exhibited discotic nematic and hexagonal 
columnar phases [351433]. The electrostatic interaction 
energy between two such dipolar ellipsoids is given by 
2 2 

U?* = E -X- [(Aia • A#0 - 3(£»<* • r a/3 ){p, jt 3 ■ f a p)\ 

a,(3=1 r «/3 

, where r a p(= r jp — r ia ) is the vector joining the two 
point dipoles fii a and fijp on the molecules i and j 
at the positions ip a = iq db and r jp = rj =b 
Then the total interaction energy between two dipolar 
molecule is given by LT° tal = U^ B + Uf ^. Here we 
have used differ ent valu es of the reduced dipole mo¬ 
ments /i* = \J /i 2 /eo<7o = 0.2, 0.4, 0.6, 0.7 and 0.9 to 
investigate the effect of dipole strength on the phase 
behavior of the system of interest. The dipole moments 
fjb* =0.2 and 0.9, for a molecular diameter of cro ~ 10A 
and an energy term eo = 5 x 10“ 15 erg corresponds to 
0.45 D and 2 D respectively. 

The Ewald summation (ES) and reaction field (RF) 
methods are two highly popular methods to handle 
the long range dipolar interaction in condensed matter 
systems. In the ES method 0 , the central cubic 
simulation cell is surrounded by its replicas in all 
directions to form an infinitely large sphere of periodic 
images. This infinitely large sphere is itself embedded 
in a continuous medium of dielectric constant e s . The 
expression for the overall dipolar interaction energy of 
the system in the ES technique is given by 
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where -F(fc) = ’ s ^{k.Hi)exp{ik.ri) and M = y^ fij. 

i i 

where ry is the position vector of the dipole moment 
Hi- M is the total dipole moment of the cubic central 
box of volume V=L 3 , erfc is the complementary error 
function and the sum on n is over lattice vectors. The 
prime in the sum over n=(n x , n y , n z ), n x ,n y ,n z integers, 
indicates that i j± j for \n\ = 0. The reciprocal space 

27r 

vectors are of the form k — —~n. The convergence of the 

1j 

interaction terms in the ES depends on the parameter 
a which is chosen such that the real space sum over 
lattice vectors can be restricted to the first term in 
the sum i.e. the term with \n\ = 0. The last term 
in the above equation for Ues, accounts for the work 
done against the depolarizing field due to the surface 
charges induced on the spherical boundary. This term 
vanishes only for e s oc (conducting / tin-foil boundary 

condition) and makes its maximum contribution when 
e s sb 1 ( vacuum boundary condition). The contribution 
to this term from the pair interactions is given by 
47r 

--——^ • If is to be noted that the term does 

(2e s + 1)L 

not decay with separation and it is clear from in¬ 
dependence that this term will favor antiparallel rather 
than parallel dipolar orientations. The choice of e s = 1 
effectively results in surrounding the large spherical 
sample of cubic replicas with a vacuum and is extreme 
for a ferroelectric medium. Therefore we can expect to 
observe a sig nificant effect upon the ferroelectric order 
found in [20] when e s is set equal to 1. The Ewald sum 
parameters used in the present work for a system of N 
= 500 dipolar ellipsoids are : ( |n maa ,| = 3, n 2 max = 27 
containing 334 k vectors , a = 15 ). 

In the reaction field method [21], the dipolar pair 
interactions are exactly evaluated upto a cutoff ra¬ 
dius Rc and the dipoles outside cutoff sphere are 
considered as a continuous medium with dielectric 
constant e s . Polarization of the continuous medium 
in response to the polarization from the total dipole 
moment M within cutoff sphere produces an electric 
field Ri that acts back on the i th dipole, located at 
the center of the cavity. The reaction field is expressed 
2(e s -l) 


as Rj = 


The size of the reaction field 


(2e s + l )R c 3J 

along any molecule is proportional to the moment of 
the cavity surrounding that molecule. The contribution 
to energy from reaction field is — which is a 

negative contribution if Hi is parallel to M. Therefore, 
the reaction field tends to polarize the system when a 


ferroelectric phase is formed, by reducing the potential 
energy of the system in favor of global ferroelectric order. 

The value of e s strongly influences the behavior of a 
ferroelectric system. If e s oc (conducting or tin-foil 

boundary condition), the contribution to energy from re¬ 
action field strongly favors a polarization. As mentioned 
before, the systems of dipolar spheres exhibited global 
ferroelectric order under the influence of the conducting 
boundary conditions in both the ES d, [.9, E-HH and 
RF methods [38], [39j. The present system of dipolar 
ellipsoids was also found to generate global ferroelectric 
order in our previous works [20|, under the influence of 
the conducting boundary conditions. If e s = 1 (vacuum 
or spherical cutoff boundary condition), Rj =0, so only 
pairwise interactions within cavity are considered. This 
boundary condition is used to simulate the behavior of a 
macroscopic sample placed in vacuum and acts against 
the formation of global polarization. If e s = 1, the 
system tends to form local structures while keeping the 
total polarization zero. The RF method with e s = 1 has 
been satisfactorily used in studies of the orientationally 
ordered systems of dipolar spheres [28], [H, [38j and model 
systems of water [4l|, (42| . 

It should be noted that in the case of ES a depolarizing 
field appears for e s = 1, which is absent for e s oc. In 
the RF method, a polarizing field appears for e s oo, 
absent when e s = 1, with the same effect. So, we expect 
to observe qualitatively similar phase behavior using 
both the methods. 

We have used three different system size : N = 1500 
(Rc = 3cr 0 ), N = 4000 (Rc = 5cr 0 ) and N = 8500 
(Rc = 7<to) when investigating the phase behaviour of 
the system using spherical cutoff approach ( e s = 1, 
Ri = 0 ) . 


III. SIMULATION DETAILS 

We have performed Monte Carlo (MC) simulation 
studies in the isothermal-isobaric ( constant NPT ) en¬ 
semble with periodic boundary conditions imposed on 
the systems of dipolar ellipsoids. We have performed a 
cooling sequence of simulation runs along an isobar at 
fixed pressure P*(= Pa^/eo) = 100. A MC simulation 
of the system of GB ellipsoids without dipoles yielded a 
discotic nematic and hexagonal columnar phases at the 
same pressure without any ferroelectric order [35|. We 
started the simulation from an well equilibrated isotropic 
liquid phase in a cubic box. We then reduced the tem¬ 
perature of the system sequentially to explore the phase 
behavior. At a given temperature, the final equilibrated 
configuration obtained from the previous higher temper¬ 
ature was used as the starting configuration. The sys¬ 
tem was subjected to long equilibrium runs at each state 
point [p*, T*]. During a MC cycle, each particle was ran- 
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FIG. 1. (color online).(a) Evolution of the average second-rank order parameter (P 2 ) against reduced temperature T* at 
constant pressure P* = 100.0 for //=0.20 and N = 1500, 4000. (b) Evolution of (P 2 ) against T* at P* — 100.0 for /i*=0.40 
and N =* 1500, 4000. (c) Evolution of (P 2 ) against T* at P* = 100.0 for /x*=0.60 and N — 1500, 4000 and 8500. (d) Evolution 
of (P 2 ) against T* at P* = 100.0 for /i*=0.90 and N — 4000 and 8500.. 


domly displaced and reoriented following metropolis cri¬ 
teria where the reorientation moves were performed us¬ 
ing Barker-Watts technique j2l|. An attempt to change 
the volume of the cubic box was also performed in each 
MC cycle. The acceptance ratio of the roto-translational 
moves and volume were adjusted to 40%. To overcome 
any possibility of locking in a metastable state, the parti¬ 
cles were also allowed to attempt up-down flip moves ex¬ 
changing particle tip with bottom with a 20 % frequency 
with respect to the roto-translational MC moves. We 
have also used an orthogonal box at some state points. 


IV. RESULTS & DISCUSSIONS 

In order to characterize different phases of the sys¬ 
tem, various order parameters were computed. The 
average orientational order of the particles was moni¬ 
tored by the second-rank orientational order parameter 
P 2 defined by the largest eigenvalue of the order tensor 
Sa0 = \(3u ia Ujp - 5 a p), where a, /3 = x,y,z 

are the indices referring to three components of the unit 
vector u along the orientation of the particles and 5 a p is 
the Kronecker delta symbol. The value of P 2 is close to 
zero in the isotropic phase and tends to 1 in the highly 
ordered phases. The global ferroelectric order was mea¬ 
sured by calculating the average polarization per par¬ 
ticle Pi defined by P\ = Yl!i=i fii-d where d is the 
director of the system. P\ is unity in a perfectly fer¬ 
roelectric phase and zero in an antiferroelectric phase 
and in the isotropic phase. P 2 is therefore the indi¬ 
cator of global orientational order and Pi distinguishes 
between ferroelectric and antiferroelectric phases. We 
have also measured the biaxial order parameter (R 2 2 ) = 
(|(1 + cos^ 0) cos 2 a cos 2 y — cos /3 sin 2a sin 2 y) as de¬ 
scribed in j43j, where a ,/?,7 are the Euler angles giv¬ 
ing the orientation of the molecular body set of axes 
w.r.t. the director set of axes. In a biaxial phase, the 
anisotropic molecules exhibit additional orientational or¬ 
der along a second macroscopic direction perpendicular 
to the primary director. It means that we can define a set 
of perpendicular macroscopic axes of preferential orienta¬ 


tion (only two need to be defined as the third is then spec¬ 
ified as perpendicular to the other two) in a biaxial phase. 
In the present system, the symmetry axes of the ellipsoids 
remain aligned in all the uniaxial and biaxial ferroelec¬ 
tric discotic phases but the molecular x axes ( axes along 
the separation between two dipoles on each ellipsoid) and 
molecular y axes ( axes perpendicular to both the molec¬ 
ular symmetry axis and molecular x axis) are significantly 
oriented only in a biaxial phase. The biaxial order pa¬ 
rameter (iil 2 ) measures the degree of ordering of the 
molecular x and y axes in a plane perpendicular to the 
primary director [43|. For (P 2 ) = 1, (^ 2 , 2 ) = 0 the sys¬ 
tem is perfectly uniaxial and for (P 2 ) = 1 , (P| 2 ) = 1 the 
system is perfectly biaxial. In order to verify the fluidity 
of the ferroelectric phases, we calculated the mean square 
displacement as follows : ( R 2 ) r = X^i[u( r ) ~ u( 0)] 2 
, where 15 (r) is the position vector of the i th particle 
after completion of r MC cycles. In the fluid phases, 
the mean square displacement steadily increases with 
increasing r indicating fluid behavior. In contrast for 
solids (P 2 ) r becomes constant as r increases. For a 
proper structural analysis of the resultant ferroelectric 
phases, we calculated important distribution functions 
as required. We have measured the radial distribution 
function g(r) = j^pr- p {5{r - rij))ij ,where the average is 
taken over all the molecular pairs. The columnar distri¬ 
bution function g c (r*) and the perpendicular distribution 
function g(r* ) for the disklike ellipsoids were calculated 
following |32j |. 


A. Results obtained using the spherical cutoff 
approach 

For /i* = 0.20, the strength of dipolar interaction 
remains quite weak to significantly influence the phase 
behavior of the GB ellipsoids. The system in this case 
exhibits a phase behavior qualitatively similar to that 
of the non-polar GB disks as reported in [35j. The 
variations of the average orientational order parameter 
(P 2 ) against the reduced temperature are presented in 
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FIG. 2. (Color online). Snapshots of the configurations generated by MC simulations at ( N=4000, /i* = 0.20) : (a) Isotropic 
phase at T* = 7, (b) Nematic phase at T* = 6 , (c) Side view of the Hexagonal Columnar phase at T* = 5, (d) The color 
plate showing different colors used in the snapshots given in this paper. The oblate ellipsoids are color coded according to their 
orientation with respect to the phase director ranging from parallel [ green (light gray)] to antiparallel [ blue (dark gray)] and 
the intermediate colors indicate intermediate orientations. All the snapshots are generated using the graphics software QMGA 

S3. 



FIG. 3. (color online). Distribution functions for ^*=0.20 (N = 4000). (a) Radial distribution function g(r*) at three different 
temperatures: T* = 6.5 (I) ,T* = 6.0 (N) , T* = 5.0 (Col); (b) Columnar distribution function g c (r jj) at three different 
temperatures: T* = 6.5, T* = 6.0 and T * = 5.0; (c) Perpendicular distribution function g(rj_) at three different temperatures: 
T* — 6.5, T* = 6.0 and T* = 5.0. The symbols I stands for isotropic, N stands for nematic phase, Col stands for the columnar 
phase. 
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(a) (b) (c) 


FIG. 4. (color online). Snapshots of the configurations generated by MC simulations at ( N=4000, /jl* = 0.40) : (a) Nematic 
phase at T* = 6.5 , (b) Top view of the striped Hexagonal Columnar phase with periodically arranged domains of nearly equal 
and opposite polarizations at T* = 6,(c) Side view of the striped Hexagonal Columnar phase at T* = 6. The particles are color 
coded according to their orientation with respect to the phase director ranging from parallel [ green (light gray)] to antiparallel 
[ blue (dark gray)]. In the columnar phase, most ellipsoids are colored either in green (light gray) or blue(dark gray) to indicate 
that most of them are oriented either parallel or antiparallel to the director. The different polarized domains can be easily 
distinguished by their colors. 



(a) 


(b) 


(c) 


FIG. 5. (color online). Distribution functions for /x*=0.40 (N = 4000).(a) Radial distribution function g(r*) at three different 
temperatures: T* = 7 (I) , T* = 6.5 (N),T* = 6.0 (SCol); (b) Columnar distribution function g c (rj|) at three different 
temperatures: T* m 7, T* = 6.5 and T* = 6; (c) Perpendicular distribution function g(r*[) at three different temperatures: 
T * = 7, T* = 6.5 and T * = 6. The symbols I stands for isotropic, N stands for nematic phase, SCol stands for the striped 
columnar phase. 



r * 

(a) 




(b) (c) 


FIG. 6. (color online). Distribution functions for /x*=0.60 (N =s 8500).(a) Radial distribution function g(r*) at three different 
temperatures: T* = 8.5 (I) , T* = 8 (SN),T* = 7.5 (SCol); (b) Columnar distribution function g c (r jf) at three different 
temperatures: T* = 8.0, T* — 7.5,T* = 6.5 (SCol); (c) Perpendicular distribution function g(r*\_) at four different temperatures: 
T* = 8.0, T* = 7.5 and T* = 6.5. The symbols I stands for isotropic, SN stands for the striped nematic phase, SCol stands for 
the striped columnar phase. 
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FIG. 7. (color online). Snapshots of the configurations generated by MC simulations for /jl* = 0.60 : (a) striped nematic phase 
with polarized domains at T* = 8 for N — 4000, (b) side view of the columnar phase with polarized domains at T* = 7.5 for 
N = 4000, (c) striped nematic phase with polarized domains at T* = 8 for N = 8500, (d) Side view of the striped Hexagonal 
Columnar phase at T* = 7 for N=8500. The oblate ellipsoids are color coded according to their orientation with respect to the 
phase director ranging from parallel [ green (light gray)] to antiparallel [blue (dark gray)]. 



FIG. 8. (color online).Snapshots of the interesting configurations generated by MC simulations for /i* — 0.90 : (a) striped 
nematic phase with polarized domains at N — 4000, T* = 11. 5,(b) striped nematic phase with polarized domains at N — 
8500, T* = 12, (c) Side view of the striped biaxial phase with domains of equal and antiparallel polarization at N = 4000, T* = 
11, (d) Top view of the striped biaxial phase at N = 4000, T* = 11 where the orientations of the dipolar separation vectors of 
each ellipsoid is indicated by red lines. The oblate ellipsoids are color coded according to their orientation with respect to the 
phase director ranging from parallel [ green (light gray)] to antiparallel [ blue (dark gray)]. 



(a) (b) (c) 


FIG. 9. (color online).Distribution functions for /i*=0.90 (N = 4000). (a) Radial distribution function g(r*) at three different 
temperatures: T* = 12 (I) , T* = 11.5 (SN),T* — 11.0 (SB); (b) Columnar distribution function g c (rjj) at three different 
temperatures: T* = 12, T * = 11.5 and T * = 11.0; (c) Perpendicular distribution function g(r±) at three different temperatures: 
T* = 12, T* — 11.5 and T* = 11.0. The symbols SN stands for the striped nematic and SB stands for the striped biaxial phase. 
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Fig. |l(a)| for two systems of different size. It can be seen 
that the behavior remain quite similar for N=1500 and 
N=4000. (P 2 ) remains close to zero at high temperatures 
identifying the isotropic phase. Upon cooling (P 2 ) jumps 
to (P 2 ) ~ 0.71 at T* = 6 indicating a transition to 
the orientationally ordered discotic nematic phase. At 
T* = 5, we observe a discotic columnar phase without 
any well defined polar arrangement of the ellipsoids. 
The snapshots of these resultant configurations are 
shown in Fig. [2] and the related distribution functions 
are presented in Fig. [3] for N=4000. The flatness in the 
radial distribution function g(r*) at T* = 6.5 reflects the 
structurelessness of the isotropic liquid in the long range. 
At T* = 6, g(r*) shows the absence of any long-range 
positional order in the nematic phase. Considerable 
structure in g(r*) for T* = 5 indicates the formation of 
the higher ordered columnar structures. The columnar 
distribution function g c (r *), which is a measure of 
positional order within a single column, is shown in 
Fig. |3(b)| The periodic nature of g c (r *) confirms the 
periodic stacking of molecules in the columnar phases. 
At T* = 5, g c (r *) decays algebraically, indicating a 
quasi longrange columnar order. The algebraic decay of 
the peak values in the columnar distribution function is 
a typical signature of columnar liquid crystal order [32| . 
For a crystal, the peak intensities remain finite at long 
range, which is clearly not the case. As the temperature 
is lowered, g c (r jj) exhibits more ordered structure. The 


small peaks for r*,r* < 0.5 in g c {r *) and in g(r*) 
describes finite probability of short ra nge f ace-to-face 
ordering in the nematic phase. Fig. 3(c) shows the 
perpendicular distribution function g(r^), which is a 
measure of translational order in the plane orthogonal 
to the orientation of the disc shaped molecules. g(r j_) at 
T* = 6.0 remains essentially flat for the nematic phase. 
At T* = 5, a flat second peak indicates hexagonal 
columnar packing. The second peak is usually broken 
for for a perfect hexagonal order. The biaxial order 
parameter (P 2 2 ) and the global ferroelectric order 
parameter (Pi) remain ^ 0 in all the above described 
phases for g* = 0.20. 


For g* = 0.40, the variations of the average orienta¬ 
tional order parameter (P 2 ) against the reduced temper¬ 
atures are shown in Fig. |1 (b)| for N=1500 and N=4000. 
The isotropic liquid condenses to a discotic nematic phase 
without any well defined polar ordering at T* =6.5 with 
(P 2 ) « 0.66. For T* < 6.0, the systems generate inter¬ 
esting columnar structures consisting of periodic array 
of ferroelectric slablike domains arranged in an antiferro- 
electric fashion. The st ripe d colu mnar phase at T* = 6 
is shown in Figs. |4(b)| and 4(c)[ It can be clearly seen 
that the system has split into a number of ferroelectric 
domains with alternating polarization in the columnar 
phase. The colors indicate two mutually anti parallel 
directions of polarization of the domains. Each domain 
consists of a number of axially polarized columns of the 
ellipsoids. The related columnar and transverse distri- 



FIG. 10. (color online). Mean-squared displacement ( R 2 ) T 
against MC cycles for the striped nematic phases at (T* = 
8,/x* = 0.60) and (T* = 11.5, g* — 0.90) for the systems of N 
— 4000 dipolar particles. 


but ion functions are shown in Figs. |5(b)| and 5(c)[ The 
formation of anti parallel polarized domains can be quali¬ 
tatively understood as follows. If a spherical ferroelectric 
sample is surrounded by vacuum, then charges induced 
on the surface of the sphere creates a depolarizing electric 
field with energy proportional to the macroscopic polar¬ 
ization. To minimize the configuration energy, domains 
are generated with opposite polarization which destroy 
the overall polarization. There is a competition between 
the energy decrease due to domain formation and energy 
increase due to domain wall formation. The domains are 
formed until these two balances each other. A similar sce¬ 
nario can be observed in case of ferroelectric crystals |44j 
and in model systems of dipolar spheres [16, 17]. Here, 
the striped phase shows a number of distinct ferroelec¬ 
tric domains with alternating polarization rather than 
different domains with different directions of polarization 
found for dipolar spheres 0. In addition the domains 
are characterized with a strong axial polarization. At 
T* = 6.5, the radial distribution function g(r*) in Fig. 
5 (a) | shows the absence of any long-range positional order 


m the nematic phase. Considerable structure in g(r*) for 
T* = 6 indicates the formation of higher ordered striped 
columnar structure. At T* = 6, # c (r*) decays rapidly, 
indicating a quasi long-range order along column axis. 
As the temperature is lowered, g c (r jj) exhibited more or¬ 
dered structure. Fig. |5(c) shows the perpendicular dis¬ 
tribution function g(r\_), which is a measure of transla¬ 
tional order in the plane orthogonal to the orientation of 
the disc shaped molecules. g(r±) at T* = 6.5 remains 
essentially flat for the nematic phase. At T* = 6.0, the 
broken second peak indicates proper hexagonal colum¬ 
nar packing. The biaxial order parameter (P 2 2 ) and the 
global ferroelectric order parameter (Pi) remain ^ 0 in 
all the above described phases for g* = 0.40. 


For g* = 0.60, the system directly transforms from 
an isotropic liquid to a striped nematic fluid at T* = 8 
for N = 4000 and N = 8500. However a smaller system 
(N = 1500) does not exhibit the striped nematic phase 
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at T* = 8. The snapshots of the st ripe d nem atic config¬ 
urations are shown in Figs. 7(a)| and 7(c) | respectively 
for N = 4000 and N = 8500. The fascinating periodic 
modulation of ferroelectric order obtained in the striped 
nematic phase is similar to that observed in the striped 
columnar phases found for /x* = 0.40. The snapshots 
in Fig. [7| clearly show that systems have generated 
slablike ferroelectric domains periodically arranged in 
an antiferroelectric fashion. Different colors have been 
used to indicate different domains of approximately 
antiparallel overall polarization. The striped columnar 
phase appears at T* = 7.5 for all the system sizes. The 
similarity in the phas e beha vior of the two larger systems 
can be seen from Fig. |l(c)| describing the variation of the 
average orientational order parameter with the reduced 
temperature. The snapshots also indicate that each 
domain is strongly polarized. 

The snapshots in Fig. [7] show that the domain size 
increases with the system size N in case of the striped 
nematic and columnar phases. For N=1500, each 
polarized slablike domain in the striped columnar phases 
consists of two rows and for N = 4000 and N=8500, 
the domains clearly consist of three and four rows 
respectively as can be seen from figs |7 (b)| and |7(d)| The 
gradual increase in the domain size clearly indicates that 
the consideration of longer range of dipolar interaction 
favors the ferroelectric order. In addition, the snapshots 
in Figd clearly demonstrate the absence of any long 
range structure in the striped nematic phases. In 
order to verify the fluidity of the nematic phase, we 
calculated the mean square displacement (R 2 ) r and 
the corresponding plot is shown in Fig. [10J The mean 
square displacement steadily increased with increasing r 
indicating highly fluid behavior at T* = 8. The related 
structural distribution functions for /x* = 0.60, N = 8500 
are shown in Fig. [6j At T * = 8, the radial distribu¬ 
tion function g(r*) in Fig. |6(a) shows the absence of 
any long-range positional order in the nematic phase. 
Considerable structure in g(r*) for T* < 7.5 indicates 
the formation of higher ordered columnar structures. 
At T* = 7.5, g c (r *) decays rapidly, indicating a quasi 
long-range order along column axis. The snapshot in 
Figs |7(b)^ md 7(d)| show the positional order along the 
columns of the columnar fluid structures. Fig. |6(c) ' 
shows the perpendicular distribution function g(r*\_) 
which at T* = 8 remains essentially flat for the nematic 
phase. At T* = 6.5, the broken second peak indicates 
proper hexagonal columnar packing. The biaxial order 
parameter (P| 2 ) an d the global ferroelectric order 
parameter (Pi) remain ^ 0 in all the above described 
phases for /x* = 0.60 

Let us now describe the phase behavior for /x* = 
0.90. The system exhibits a striped biaxial phase hav¬ 
ing slablike ferroelectric domains in addition to the uni¬ 
axial striped nematic phase. The phase biaxiality orig¬ 
inates exclusively from the strong dipolar interactions. 


No columnar phase is obtained at this dipole strength. 
The variations of the orientational order parameter (P 2 ) 
against the reduced temperatures for different system 
sizes are shown in Fig. |l(d)[ It can be seen from 
Figs. |l(c)| and |l(d)| that the temperature range of the 
striped nematic phase gradually increases with the dipole 
strength for a fixed system size which is similar to the be¬ 
havior of the system when studied under the influence of 
the conducting boundary conditions in our previous in¬ 
vestigation [40| where the temperature range of the single 
domain ferroelectric nematic phase gradually increased 
with /x*. Here, a system size effect is also observed as 
the temperature range of the striped nematic phase is 
found to increase with N as indicated by the plots in 
|l(d)| The largest system consisting of N = 8500 dipolar 
ellipsoids exhibits a striped nematic phase at T* = 12 
with (P 2 ) ~ 0.5. The ordering of ferroelectric domains 
as shown in Figs. 8 (a) and | 8 (b)~| is quite similar to that 
observed for the striped nematic and striped columnar 
phases obtained for lower dipole strengths. Fig. |9(a) 
shows the radial distribution function g(r*) at different 
temperatures. The flatness in g(r*) at T* = 11.5 reflects 
the structurelessness of the striped nematic liquid in the 
long range. At T* = 11.5, # c (r*) and g±(r* L ) also show 
no sign of any long range positional order. We have also 
calculated the mean square displacement ( R 2 ) r for the 
striped nematic phase at T* = 11.5. The mean square 
displacement steadily increased with increasing r indi¬ 
cating fluid behavior as shown in Fig. [TOj Considerable 
structure in g(r*) at T* = 11 indicates the formation 
of a more ordered biaxial phase. We have measured the 
conventional biaxial order parameter (# 2 , 2 )? th e va lne 
which is ~ 0.95 as measured in the biaxial phases. How¬ 
ever, (i ?2 2 ) remain ^ 0 in the striped nematic phases 
reported here. The columnar and tran sverse distribu¬ 
tion functions shown in Figs. |9(b)| and 9(c)[ are quite 
different from those obtained for the columnar phases. 
From the snapshots of the biaxial phases, we can under¬ 
stand the related structure. From Fig. | 8 (d)| it can be 
seen that the dipoles are arranged such that the dipo¬ 
lar separation vectors pointing from one dipole on a disc 
to the other get oriented mostly in same direction. This 
ordering can be understood as the effect of strong pair in¬ 
teraction possible in this condition between two terminal 
dipoles of neighboring ellipsoids if two parallel dipoles are 
placed above each other. The snap shot of such interca¬ 
lated arrangement is shown in Fig. | 8 (c)[ The structural 
distribution functions for the biaxial phase are shown in 
Fig. |9(a)p(c)| At T* = 11 , g c (r*) indicates the interca¬ 
lated arrangement of the ellipsoids. The shorter peaks 
generate due to the presence of the pairs of side-by-side 
ellipsoids covering single ellipsoids. It should be noted 
that strong dipoles stabilize striped nematic phases with 
relatively weak orientational order. The global ferroelec¬ 
tric order parameter (Pi) remained approximately equal 
to zero indicating zero net polarization for all the phases 
observed for /x* = 0.90. 
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B. Results obtained using the Ewald summation 
method 


We have also studied the phase behavior using the con¬ 
ventional but more time consuming Ewald Summation 
method with e s = 1. We have chosen a smaller sys¬ 
tem consisting of N = 500 dipolar ellipsoids to study the 
depolarizing effects using ES approach and investigated 
the phase behaviors for two different dipole strengths 
/x* = 0.70 and /x* = 0.90 for which we observe the fer¬ 
roelectric nematic, ferroelectric columnar and ferroelec¬ 
tric biaxial order within periodic slablike ferroelectric do¬ 
mains. 


For /x* = 0.70, the striped nematic phase with 
ferroelectric slablike domains is obtained at T* = 9 
upon cooling the isotropic phase at the fixed pressure. 
The snapshot in Fig. |12(a) clearly shows the pres¬ 
ence of ferroelectric nematic domains with mutually 
opposite polarization at T* = 8.5. Further lowering of 
temperature gives rise to a columnar liquid crystal, at 
T* = 8, with ferroelectric domains consisting of three 
rows of columns. The plots of the structural distribution 
functions are presented in Fig. [13j The flatness in the 
radial distribution function g(r*) at T* = 8.5 confirms 
the absence of any long-range positional order in the 
striped nematic phase shown in Fig. |12(a)| Considerable 
structure in g(r*) for T* < 8 indicates the formation of 
the higher ordered columnar structures. The columnar 
distribution function g c (r *) is shown in Fig. |13(b)[ At 
T* = 8, g c (r *) decays algebraically, indicating a quasi 
longrange columnar order as shown in the snapshot in 
Fig, 12(b) The small peaks for r*,r* < 0.5 in g c (r*) 
and in g(r*) indicates finite probability of short range 
face-to-face ordering in the nematic configurations. Fig. 


13(c) shows the perpendicular distribution function 


at T* = 8.5 remains essentially flat for the nematic 
phase. At T* = 8, a broken second peak indicates 
hexagonal columnar packing of the columns. The above 
characteristics are qualitatively similar to that of the (N 
= 4000, /x* = 0.60) system investigated using the SC 
technique. The biaxial order parameter (# 2 , 2 ) an d the 
global ferroelectric order parameter (Pi) remain « 0 in 
all the above described phases for /x* = 0.70. 


For /x* = 0.90, the striped nematic phase having 
ferroelectric slablike domains is obtained at T* = 11 
upon cooling the isotropic phase. The variations of 
the average orientational order parameter (P 2 ) against 


the reduced temperatures are presented in Fig. 12(e) 


and the snapshots of the interesting configurations 
obtai ned using t he Ewald sum method are given in 
Figs. |12(c)p2(d) The presence of ferroelectric nematic 
domains with antiparallel polarizations at the above 
temperature can be seen from the snapshots. Further 
lowering of temperature generates a highly ordered 
biaxial phase, at T* = 10.5, with ferroelectric domains 
of width « 3 <to. The strong dipolar interactions in this 
case stabilizes a biaxial phase instead of a columnar 


phase. The characteristics are again found qualitatively 
similar to that exhibited by (N = 4000,/x* = 0.90) 
system studied using the spherical cutoff approach. 
The plots of the structural distribution functions are 
given in Fig. [14l The flatness in the radial distribution 
function g(r*) at T* = 11 confirms the absence of any 
long-range positional order in the striped nematic phase. 
Figs. |14(b)| and |14(c)| show that the columnar and 
perpendicular distribution functions remains essentially 
flat for the nematic phase at T* = 11. Considerable 
structures in g c (r*) and g(r*_ j_) at T* = 10.5 indicate 
the formation of the higher ordered biaxial phase. We 
have measured the biaxial order parameter (# 2 , 2)5 the 
value which is ~ 0.95 as measured in the striped biaxial 
phase at T* = 10.5. However, (#2 2 ) remain ^ 0 in the 
striped nematic phase. The columnar and transverse 
distribution functions at T* = 10.5 are found quite 
different from those obtained for the columnar phases. 
From the snapshot of the biaxial phase in Fig. |12(d)[ 
we can understand that the structure is quite similar 
to the biaxial phase obtained using the spherical cutoff 
approach. It can be seen that the dipoles are arranged 
such that the dipolar separation vectors pointing from 
one dipole on a disc to the other get oriented mostly in 
same direction. This ordering can be understood as the 
effect of strong pair interaction possible in this condition 
between two terminal dipoles of neighboring ellipsoids 
if two parallel dipoles are placed above each other. At 
T* = 10.5, g c (r*) indicates the intercalated arrangement 
of the ellipsoids. 

The formation of a biaxial phase via dipolar inter¬ 
actions supports earlier simulation works on dipolar 
anisotropic particles which showed that dipolar in¬ 
teractions can introduce a biaxial orientational order 
even in the absence of a shape biaxiality of the con¬ 
stituent molecules [35|, HI). The global ferroelectric 
order parameter (Pi) remained approximately equal to 
zero indicating zero net polarization for all the phases 
observed for /x* = 0.90. We have observed that the 
system developed a small cavity at this dipole strength 
which might be a result of the fixed shape constraint of 
the cubic box used by us as reported earlier in studies of 
discotic liquid crystals in (37]. 

Therefore, the simulations of smaller systems per¬ 
formed using the ES technique, generated results qual¬ 
itatively similar to that obtained from the simulations of 
larger systems using the spherical cutoff technique. Usu¬ 
ally it is considered that the macroscopic properties of 
a system of dipolar particles converge more rapidly with 
systemsize in case of ES method (e s = 00 ) than in case 
of RF method. Here, the results also indicate that the 
macroscopic properties of a system of dipolar particles ex¬ 
hibiting long range ferroelectric order also converge more 
rapidly in case of the ES method (e s = 1) than in case 
of the simpler spherical cutoff method. A study on the 
system size effects using ES technique, although of un- 
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FIG. 11. (color online). The variation of the pair interaction 
energy Ut otai of two polar ellipsoids having parallel dipolar 
separation vectors as one polar ellipsoid is slid over another 
along dipolar separation vector with their dipoles oriented 
along their symmetry axes and keeping the component of the 
inter particle vector parallel to the symmetry axes of the par¬ 
ticles at a fixed value of a e . r± is the component of the inter 
particle vector perpendicular to the symmetry axes of the par¬ 
ticles. The different curves are for different combinations of 
d* and /jl* as described inside the figure. 


deniable interest, has not been attempted here because 
of the substantial computing cost required. However, we 
think that the domain size would increase with system 
size in the ES method similar to behavior observed in 
spherical cutoff approach. A similar view was also dis¬ 
cussed for ferroelectric systems of dipolar spheres in [Tfj. 

Let us now discuss the behavior of the pair interac¬ 
tions. The pair interaction energy between two dipolar 
ellipsoids is plotted in fig. [Tl] as a function of the 
component of the pair separation vector r perpendicular 
to the particle symmetry axes as one particle is slid 
over the other maintaining a fixed separation (r M = cr e ) 
along the symmetry axis and having all dipoles oriented 
in same direction. Note that when r\_ is zero one 
particle completely covers the other. Curves are shown 
for d* = 0 and d* = 0.50. Here we have discussed the 
case where the ellipsoid is sliding along the direction 
of molecular x axis, i.e. along the dipolar separation 
vector d*(= d*x ) of the ellipsoids. This particular case 
is important for the biaxial phase behavior found for 
/i*=0.90. It is evident from Fig. [TT| that the dipolar 
separation d* has a dramatic role on the pair potential. 
For d* = 0, the strong head-to-tail dipolar interaction 
results in a sharp minima at r* L = 0. It is the sharpness 
of this potential which leads to strong stabilization 
of the columnar phase for central dipolar disk shaped 
particles studied in [9. [46|. In the columnar phases of 
central dipolar disks, the individual polarized columns 
do not interact strongly enough and entropy ensures that 
equal numbers will be polarized in opposite directions 
P, |46j. In the present model, the pair interaction looses 
its sharp focus when the dipoles are more separated 


and generates a weaker and slowly varying attraction 
which helps in generating the ferroelectric nematic order 
for d*= 0.50. It can be seen that for separation d* = 
0.50, a second well is constructed with its minimum 
at r^=0.50. The second well simply results from the 
interaction between two parallel dipoles when they 
are positioned above each other. The strengthening 
of this second well favors strong stabilization of the 
ferroelectric biaxial phase at higher dipole strengths. 
In the uniaxial striped nematic and striped columnar 
phases, the dipolar separation vectors of two neighboring 
ellipsoids are usually randomly oriented with respect to 
each other. In such cases, the pair interaction do not 
generate a second well as found for parallel orientation 
of the dipolar separation vectors of the two ellipsoids. 
Here, the formation of the ferroelectric columnar do¬ 
mains for /j,* = 0.40,0.60 and 0.70 indicate that the 
column-column interactions are stronger than that in 
case of central dipolar disks. 


V. CONCLUSIONS 


In this work, we have systematically investigated the 
influences of the depolarizing boundary conditions upon 
the existence of a class of novel ferroelectric phases of 
dipolar origin in the systems of disklike anisotropic par¬ 
ticles using both the spherical cutoff and the Ewald Sum¬ 
mation techniques. The systems resulted in the forma¬ 
tion of periodic slablike ferroelectric fluid domains ar¬ 
ranged in an antiferroelectric fashion. We found exis¬ 
tence of the long range ferroelectric nematic, ferroelec¬ 
tric columnar and ferroelectric biaxial orderings within 
the domain regions at different state points and dipole 
strengths. It is found that the minimum value of the 
dipole strength required to generate a ferroelectric or¬ 
der is fi* = 0.40. The periodicity and the size of the 
slablike domains, for a fixed systemsize, remained un¬ 
changed with respect to a variation of temperature or 
dipole strength studied here. It is explicitly shown that 
the size of the ferroelectric domains grow with the sys¬ 
tem size as an effect of considering increasing range of 
dipolar interactions. These results provide great support 
to the novel understanding that the dipolar interactions 
are indeed sufficient to produce a long range ferroelec¬ 
tric order in discotic fluids. Our results also supports the 
view expressed in [16] that for sufficiently large samples 
the local orientational order in the ferroelectric domains 
would be quite similar to that obtained with e s = oo. Our 
study should be helpful in analyzing and understanding 
the structural properties of novel anisotropic fluids where 
dipolar forces would play a dominant role in generating 
global ferroelectric order. 
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FIG. 12. (Color online). Snapshots of the configurations generated by MC simulations at N = 500 using Ewald Summation 
method : (a) striped nematic phase at /jl* — 0.7, T* = 8.5, (b) Striped columnar phase at /jl* — 0.7, T* = 7.5 , (c) Striped 
nematic phase at fi* = 0.9, T* = 11, (d) Striped biaxial phase at /jl* = 0.9, T* = 10.5. The oblate ellipsoids are color coded 
according to their orientation with respect to the phase director ranging from parallel [ green (light gray)] to antiparallel [ blue 
(dark gray)] and the intermediate colors indicate intermediate orientations. 
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FIG. 13. (color online). Distribution functions for /jl* = 0.70 (N = 500). (a) Radial distribution function g(r*) at three different 
temperatures: T* = 8.5 (SN) ,T* = 8.0 (SN) and T* = 7.5 (SCol); (b) Columnar distribution function # c (rjj) at three different 
temperatures: T* = 8.5,T* = 8.0 and T* =s 7.5; (c) Perpendicular distribution function g(r±) at three different temperatures: 
T* = 8.5, T* = 8.0 and T* = 7.5. The symbols SN stands for striped nematic phase, SCol stands for the striped columnar 
phase. 
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FIG. 14. (color online). Distribution functions for /i*=0.90 (N = 500). (a) Radial distribution function g(r*) at two different 
temperatures: T* = 11 (SN) and T* = 10.5 (SB); (b) Columnar distribution function g c (r jj) at two different temperatures: 
T* m il and T* = 10.5; (c) Perpendicular distribution function g(r±) at two different temperatures: T* = 11 and T* — 10.5. 
The symbols SN stands for striped nematic phase, SB stands for the striped biaxial phase. 
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